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ABSTRACT 



Prediction of demand and setting reorder levels by ex- 
ponential smoothing and mean absolute deviation is studied 
for an inventory system. The relation between the stand- 
ard deviation and mean absolute deviation of demand dis- 
tribution is discussed. Asymptotic first and second 
moments of the exponentially-smoothed forecast errors (MAD) 
are derived. Results obtained from simulation of several 
normal systems indicated that the smoothing technique is 
inferior to the classical maximum likelihood estimation 



method. 
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TABLE OF SYMBOLS AND ABBREVIATIONS 



th 

Demand during the i time period 
Oi Exponential smoothing constant 

til 

S^(x) Forecasted demand during the i^ period 

p Correlation coefficient 

0 Risk level 

A Mean absolute deviation of a probability 

distribution 

fl Mean demand 

V( ) Variance of the argument 

E( ) Expected value of the argument 
til 

e^ t^ period forecast error 

CT Estimate of the standard deviation of the 

demand by maximum likelihood estimates 
cr Estimate of the standard deviation of the 

demand distribution by mean absolute 
deviation 
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I . INTRODUCTION 

In recent years the technique of forecasting time 
series by exponential smoothing has received a considerable 
amount of attention in the literature. As a forecasting 
technique in inventory theory, it is particularly suited to 
the situation where demands are time dependent and corre- 
lated for much of its intuitive appeal rests on weighting 
the past in such a way that the most distant past receives 
less and less weight. When demands are not so related, 
smoothing may or may not be a suitable forecasting tech- 
nique. In particular, when demands are independent and 
from a common distribution such as for the so-called fast 
moving items, more classical methods like maximum likeli- 
hood may be preferable because of the optimal properties 
that are often demonstrable. Even in such circumstances, 
however, exponential smoothing is still a candidate for a 
forecasting scheme and it may very well be that, while not 
optimal , the difference in such criteria as stockout risk 
for example, may be negligible compared to optimal proce- 
dures from a practical point of view. In that case, fore- 
casting rules need not be changed for items differing in 
demand pattern and a certain uniformity may thus be 
achieved. 

The main purpose of this thesis is to investigate a 
model of an inventory system used by the Naval Supply 
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System Command (NAV SUP) following Brown [Ref. l]. The 
model, as described by Ref. 1 and Ref. 5, is briefly as 
follows ; 



Single exponential smoothing is an operator given in 
Ref. 1 as 



a is the smoothing constant with 0 < a < 1 

/9 = 1 - a 

Intuitively speaking the smoothing operator is a 
method of forecasting by weighting past demands in such a 
way that the most recent demand is weighted by OL and the 
remaining demands are given less and less weight as time 
progresses . 

The standard deviation, CT, of the demand distribution, 
which is needed to set the reorder levels, is estimated by 
the following three statistics; 

e^ = - S^_^(x), t*"^ period forecast error. 



t-1 




( 1 - 1 ) 



k=0 



til 

where x^ is the total demand of the i period 



t-1 




( 1 - 2 ) 



k=0 
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estimate of the mean absolute deviation of forecast errors 



and 




estimate of the standard deviation of the demand distribu- 
tion. 

This estimate of o, as explained in Ref. 1, is based 
on the mean absolute deviation (MAD) , A, of the forecast 
errors. Reorder levels, based on the predetermined risk 
level j|) (defined to be the probability of running out of 
stock), are found as follows: 

0 = P(X > |i + kCT) 

Here, 4 is the mean of X, and k is a constant which 
can be determined from the probability tables for a given 
risk level, ij). The reorder level is then given by 

Reorder Level = 4 + kCT (1-3) 

where k is the given above. 

In the analysis this model was slightly modified be- 
cause of the following considerations: 

1. X , the demand asked from the inventory system 
o 

before it is placed in operation, is zero. Although one 

may predict a number for a fictitious inventory, there will 

be no observation for X . Hence e is zero. 

o o 
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2. 0^ and can be assumed to be approximately equal 

to zero for values of N > 30, if 0.01 s a s 0.99. In view 
of this, the modified form of (1-1) and (1-2) used in the 
analysis was 

t-1 

Sj.(x) = a Y 

k=0 

t-1 

A = a Y ^ ^ ®t-k' 

k=0 

In this model, not only is demand forecasted by expo- 
nential smoothing but demand variability is estimated by 
means of exponentially smoothed absolute deviations. Such 
estimates are in turn used to set safety levels for inven- 
tory order rules. 

The model described above has been subjected to many 
criticisms. Zehna [Ref. v] has indicated that the estimate 
of standard deviation might be very poor depending upon the 
underlying demand distribution. Based on the results ob- 
tained by Zehna [Ref. s], the analysis in this paper takes 
on its starting point on investigation of the ratio A/a, 
a key element in the estimation of a. It is quite evident 
that any weakness in the estimation procedure for a affects 
the reorder levels. An underestimate imposes higher risk 
and an overestimate causes unnecessary expenditures. In 
particular, the resulting consequences could be catastroph- 
ic in a military inventory system. 
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It is quite evident that the probability distributions 
of certain key statistics would be desirable in any exten- 
tion of the existing analysis. Although some of the dis- 
tributions required have been derived, the distribution of 
A , a most important statistic, could not be derived ex- 
plicitly. However, the first and second asymptotic moments 
of A were obtained, and it is shown that the asymptotic 
variance of A is finite for the case where the demand has 
a normal distribution. Therefore, the next step of the 

analysis used simulation techniques to gain an insight into 

❖ 

the behavior of the distribution of A , and to observe the 
effects of a finite variance. Reorder levels and their 
dispersions along with the risk levels attained are analyz- 
ed. The resulting errors due to the variability of safety 
levels is examined in detail and the consequences are 
evaluated from several points of view when demands are 
normal. Moreover the results are compared to maximum like- 
lihood methods establishing that a considerable difference 
exists in the two methods. For items satisfying the demand 
assiimptions considered here, then, it is established quan- 
titatively that exponential smoothing is definitely inferior 
to more classical methods. 
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II. RATIO OF A/g 



The mean absolute deviation of a probability distri- 
bution is given by 

A = E(| X - 4 I ) . 

For a normal distribution the ratio A/ct is the con- 
stant A/2/ff’ = 0.8. Brown [Ref. l] in his work on smoothing 
and prediction states, 

"Hence, the mean absolute deviation is proportional to 
the standard deviation, and the ratio may depend on the 
form of the distribution p(t), but dependence is slight." 

Now, as indicated by Zehna in [Ref. 7 ], this is not 
always the case. For example, consider a Poisson distri- 
bution with parameter X. A/c is then a function of X and 
it may be much smaller than 0.8 for small values of X 
[Ref. 8]. However, there are some distributions where the 
approximation 0.8 for the ratio A/a is a good one under 
certain conditions. The quote above is a rather general 
statement which may be justified as follows; 

Let X be a continuous random variable with probability 
density function (p.d.f.) f(x) and let u and _1 be the upper 
and lower limits, respectively, of the random variable X on 
the real line. Then 

u 

A = E(lx-fi|) = f|x-u| f(x)dx 

•’a 
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or 



p|i nU pU pU 

xf(x)dx + (i f(x)dx + xf(x)dx - fj. f (x) dx. 
^ Z £ fl \ 



Hence, 



A = 2 

Here, 



p4 

F(x)dx 



F(x) is the cumulative distribution of X. 



Then 



A 

<7 



g 

a 



F(x)dx 



Let X = ory + )i. Then, 



A 

a 



2 



4 

f 

F(cry +(i) dy 
£-U 

a 



Now suppose that F(x) is any distribution function such 
that 



F(x) 

where ^(z) 
function. 



is the cumulative standard normal distribution 
Then for large values of ju/cT 



F(CTy + m) S $(y) 
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Hence 





dz dy 



Integrating by parts, the ratio then becomes A/o’ = / 2/ff'. 
There are some distributions for which the above approxima- 
tion holds; for example erlang and chi-square. But for 
other distributions such a generalization may not apply. 

For instance, let X have the following mass function 



X = 0 
X = a 
X = 2a 

Then 

A = 0.2a 

O^ = 1.2a 

and 

A/o = .182. 



p(x) 



^0.1 


if 


00 

• 

o 


if 


Vo.i 


if 



The error made by assuming A/o = 0.8 is 0.618, or 11 % , 
which is not negligible. 



III. PROBABILITY DISTRIBUTION OF A 
AND OTHER ASSOCIATED STATISTICS 

Of central interest in the use of MAD is the ability 

to estimate the demand variability. The variance of the 

forecast error is always a function of the variance of de- 
2 

mand, CT . Therefore the MAD of the forecast errors is also 
a function of cr. Since MAD is a linear function of the 
standard deviation [Ref. l] and if MAD of the forecast 
errors is estimated, then cr can be estimated as described 
in [Ref. l]. But in this estimation there are two points 
to be considered. 

1. Invariance Principle ; 

CT and A are related to each other by a constant, 

i . e. , 

CT = cA . 

But this relation is in the population, and it does not 
mean that the estimator of CT and A will enjoy the same 
relationship necessarily. If these estimates mentioned 
above were MLE, then the invariance principle does hold. 
However, the definition of A indicates that this is not 
the case. 

2 . The Constant of Proportionality of A/ct; 

If the underlying demand distribution is not 
normal, then the assumption that the proportionality con- 
stant is approximately 0.8 needs to be validated. For this 
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case the current model can be applied by a slight modifica- 



tion on or, that is, replacing vff/2 by the reciprocal of 
the appropriate proportionality constant. 

In this section we assume that demand is Normal (a,CT ), 

where a is a constant. Derivation of the probability dis- 

K .... 

tribution of A requires probability distribution for the 

following statistics: 

A. EXPONENTIAL SMOOTHING OF PAST OBSERVATIONS 

(Estimate of the Expected Value of Demand) 

t-1 

s^(x) - a 

k=0 



Since all X's are mutually independent, the 

distribution of S^.(x) for large values of t is known to be 

2 01 

a normal distribution with mean a and variance CT (1 + ) . 

B . FORECAST ERRORS 

t-i-2 

e,. . = X,. . - a ) ^ ^ i < t 

t-i t-i t-i-k-1 

k=0 



e^ ^ has a normal distribution with the follow- 
ing mean and variance: 

li = 0 for large values of t; 



V(e^_i) = 



(1 + 



a 



1 +^ 



for large values of t. 
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C. JOINT DISTRIBUTION OF FORECAST ERRORS 

Let U and V be t-dimensional vectors defined 



as 



®t 




1 

cj 

rr 

I 


1 — 1 
1 



(D 


u = 


^t-i 

• 

• 

• 

• 

• 

• 


• 

• 




^2 


®1 


i 





where U^, U^. are constants. Let Y be defined as 

T 

Y = U V 



or 



Y 



t 
< — » 

) U.e. 

L, 1 1 

i=l 



Substituting e^ in terms of the x's. 



t i=2 

i=l k=0 



or 



i=l 



' I “ I i ^i-k-i 

i=l k=0 
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Hence 



t-1 i-1 




i=0 k=0 



Since the coefficients of the X's are constants and the X's 
are mutually independent normal variables, it follows that 
Y has a normal distribution. Hence by the definition of a 
multivariate normal distribution [Ref. 2], the forecast 
errors have a multivariate normal distribution. The corre- 
lation coefficient between e. . and e. where i < j, can 
be computed as follows. 



Let 




where i < j < t. 



V(Y) = V(e^ .) + (V(e. .) + 2Cov(e, . , e^ . 

' t-x t-j t-i t-j 



or 




Rewriting in terms of the X's, 



Y = X 



t-i 



t-i-1 



- a0X 



t-i-2 



. . - aj9^“^"^X 



t-j 






t-j-1 



. . . + X 



t-j 



aX 



t-j-1 



- a$X 



t-j-2 
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Rearranging terms on the right-hand side 



j-i-2 

Y = x^_. - a I - 

k=0 

t--j-l 

k=0 



Since the terms in the above expression are mutually inde- 
pendent, the variance of Y for large values of t can be 
found as follows: 



V(Y) = (l + OL^ ^ ^ + (1-ajSP ^ ^ ^ ^ 

1 -^"^ 1 - $ 



Substituting V(Y), V(e^_j) into (3-1), 



.2 , 



Cove(e,.i,e^..) = f (l . ^ 2 , ^ 



-T^) 

Hence, the asymptotic correlation coefficient between 

and e, . becomes 
t-D 






(3-2) 
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Now consider the following function. 



p' = ^ 

^ 4 



(i + (1+^^)"^) + (1-0^^^ 1 



1 +^ 



Here y is a continuous variable. Then 



dy 4 









2.2y 



l+;3/^ 2oc$^ T ^ , 2a^$ ^ ^ a'\ 






Rearranging the terms. 



M', ,y ^ , I . 4 . ii±|l« , ) 



or 



M’ _ logjS (i+^)g ~A _ log ^ 

dy 2 V j3 y " 2^ 

Now log ^ is negative; therefore, the first derivative of 
p', with respect to y is positive. Hence p' , is a monotonic 
increasing function of y. Then positive integral values of 
y also lie on this same continuous curve. Actually, p co- 
incides with p' , whenever y is a positive integer. Hence p 
is a monotonic increasing function of j-i. The supremum of 
p is found by letting j-i “•*’ in (3-2) ; 

Sup p = ^ (l + (1+1) + l)- 1 

= 0 
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The minimum value of p is found by letting j-i = 1: 



Pmin = ^ IT? + <!-“)") - 1 

= ^ (l + (l+/9)a + ]3^^' - 1 

= i±^ - 1 
2 



The minimum value of $ may be zero; 
D. DISTRIBUTION OF | 

Rewriting | 




hence, p . 

min 



-0.5 . 



The term 



if?' 



has, asymptotically, the chi-square distribution with one 
degree of freedom [Ref. 3]. 

The term 



^ a' 



t -i 

( 1 + 



a 



l+$ 
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has, asymptotically, the chi distribution. Hence the mean 



value and variance of is given by 



E( 





(3-3) 



V( 




n-2 

n 



(1 + 



a 



l+$ 



(3-4) 



An attempt to derive the probability distribution of d was 
made utilizing the change-of-variable technique, and by 
defining the following transformation: 



t-1 




k=0 



Y 



1 



e 



t-1 



^t-l= ^1 



The inverse transformation becomes 



e = ± Y 
1 t-1 



^2 = ^ ^ t -2 



t-1 

-t = 

k=l 

This transformation is obviously not one-to-one; there 
are 2^ sets in the domain of the transformation that are 
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transformed into a single set in the image. For all cases 
the absolute value of the Jacobian of transformation is 
one. The joint density of A , then becomes 

^ t-1 

f(Yl» ^ 2. 1 

i=l k=l 

where g(x-, x~, ..., x ) is the p.d.f. of a multivariate 

12 n 

normal distribution, and for every i one combination of the 

/\ 

signs is valid. The distribution of A then naturally re- 
quires the (t-1) -fold integration of the joint density of 
A , Y^,..., . Although every parameter in the joint 

density is known, the labor required to arrive at an ex- 
plicitly-defined probability distribution makes the task 
very difficult. Therefore, the analysis concerning the 
dispersion of a was not completed because of the intracti- 
bility. 
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IV. MOMENTS OF A 



The next step of the analysis, then, was a simulation 
of the inventory system in order to gain an insight concern- 
ing the distribution of A . It was seen that the reorder 
levels were fluctuating around the theoretical reorder lev- 
els, even after 1000 observations. This was related to the 
asymptotically finite variance of A as indicated by [Ref.?]. 
Although computer simulations implied variance of A does 

A 

not vanish, the moments of A were derived in order that an 
analysis of the simulation results might be made. The mo- 

A 

ments of A were derived for large values of t and are, 
therefore, asymptotic. 

A 

A. MEAN VALUE OF A FOR LARGE VALUES OF t. 
t-1 

^ “ I I ®t-k I ^ ° 

k=0 

t-1 

E(A) = a ^ E( I e^._^ I ) 
k=0 




A 

B. VARIANCE OF A 



t-1 t-2 t-1 

V(A) = a^ Y e^._J )+2a^ ^ ^ ^^'^^Cov(I e^._J ,| ) 

k=0 i=0 j=i+l 

(4-1) 
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Cov(|e|._.|, |e|._.|) = E(| e^_.e^_.|) - E(| e^_^| )E(| e^_^| ) 

(4-2) 

Thus, substituting (3-3), (3-4), ancJ (4-2) into (4-1), 



‘n-2 _2 2 a. 

V(A) = -J- CT 3T5 aJ3 - 



M. 



U+P) n 



a2 + 



t-2 t-1 



2a 






(4-3) 



i=0 j=i+l 



This last expression for V(A) requires the computation of 
E(| j I ) • following derivations concern 

The joint distribution of and indicated 

in Section III, is bivariate normal. The parameters of 
this distribution are 



“i = = ° 



^1 = ” 



^2 = = 0 



/n ^ 

'^2 + IT 



1+3 



and asymptotic p is given by (3-2) . For notational con- 
venience, let 



X = e,_. 



Y = e. . 
t-3 
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Then the joint p.d.f. is given [Ref. 4] by 



f(x,y) = 



_ a 
2 



/ 2 
21T a a vl-P 

X y 



- 00 < X < 

— 00 < y < OO 



where 



5 ^—2 - 



1 ^ X 



1-p a 



X 



X y CT 



In that case 



CX> 00 



E(|xy1 ) = I xyl f(x,y)dy dx 



_00 _00 



CO Q 00 

= J I xl (^-J y f(x,y)dy + J Y f(x,y)dy^ dx 



_oo _00 



o o 



_oo _oo 



■» p p 

xy f(x,y)dydx - xy f(x,y)dydx - 



-00 o 



»00 



J J xy f(x^y)dydx + J J xy f(x,y)dydx 



o 



o o 



However, by the definition of correlation coefficient. 



o o 



E(XY) 



= xy f(x,y)dydx + xy f (x,y) dy dx + 



_oo '' _00 



-C» Q 



/xO 

J 

o 



00 00 



J xy f(x,y)dydx + J j x y f(x,y)dydx 



o o 



a a p 
X y 
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Thus 



E(1 xy| ) + E(XY) 




00 00 

f(x,y)dydx + J J xy f (x, y) dydx^ 
o o 



Making a change of variable in the first integral, let 



z = -x 
u = -y 



then 



o 


O 

ft 1 


oo oo 

r* p 


I 


xyf(x,y)dydx = 


zu f(-z,-u)dudz 


^ _oo 


— 0® 


o o 


Since 


f (x,y) = f (-x,-y) , 


, it follows that 




1 

xyf(x,y)dydx = 


00 c» 

zu f(z, u)dudz 


— 00 


^ -<x> 


o o 



Thus 



r r 

E(|xy1) + E(XY) = 4 J xy f(x,y)dydx 



(4-4) 



The right-hand side of (4-4), without the coefficient 4, is 



X 



(v-b) ^ 



J xy f(x,y)dydx = 



2CT 



X e 



O ^x^/ ^ ^ 



ye_ 



2fT^^(l-p^) 



o ^/2’^ a^yi-p^' 



dy dx 



(4-5) 

where b is the conditional mean of Y given X, and is given 

by [Ref. 4], ^ 

b = ^ PX 
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The conditional standard deviation of Y given X is given by 




Let 



y - b 



= z , then 



_ (v-h) 

2 2 
2a 2c 

xe ^ r ye 



o o 



dydx = 



X 

00 cr 00 

X 



2ie f (c2+^)e dz dx 



-b 

c 



Integrating the inner term first; 



(4-6) 



cz+b 



-b 72 IT 

c 



<iz = 



2c' 



a/27T* 



+ b (1 - 






where 4>(x) is the cumulative normal distribution. Sub- 
stituting this expression into (4-6) , 



X 



2a 



xe 



X 



■^o a 72 tt 

X 



X 



ye_ 



I 2 

_ (y-^ ) 

2c2 



X 



o C 72 it 



dydx = 



X 



2 2 

CO 2a ^ 2c 

X 

cx e 

o a X 2 tt 

2 



dx + 



2a 



2a 



bx e 



X 



o a 72TT 

X 



dx - 



b ^>(- S )e 
1 c 

■ O a 72 tt 

X 



dx 



(4-7) 
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Denote the terms on the right-hand side as 1, 2, and 3, 
respectively, and compute each one separately as follows: 



1. Substituting values of c and b into 1, 



«2 2 
P X 



X 



cx e 



^^x 2(1 -P^)(Tj,2 






2 aJ{l-p^ 



2WCT 



dx = 






x e 



2 fra 



X 



dx 






o_ cr (1-p ) 

21T X 



2 — 

a^a^{l-p^)2 

2tr 



2. Substituting the value of b into 2, 
2 



2 

per X e ^ cr p 

dx = ^ 

o O’ 

X 



o o 

X 



2o 

2 ^ 
X e 



dx 



a (7 p 
_ V X 

2 



3. Substituting the values of c and b into 3, 

2 



^ 2a 

“ ^ ' X a 

4> r_ J X e, _y 

o a.Vl-p^' a„y2ff 



X p = 



X 



X 



_2S 2 

1 __2 _ X I e 



a 7217' o 

X 



J e - J 



y2ir 



— dt dx 



where 



P = 



Px 



a^v/l-p 






(4-8) 
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Let. x/a^ = y in (4-8) ; then (4-8) becomes 



a a p 

X 



00 ^ 

2 /2 
Y e 
o 



-r 






dt dy 



where 



r = 



py 



Jl-P^ 



' f f 



Integrating by parts. 



2 

_ iL- 

du = y e , dy 



V = y 



-r - 






dt 



u = - e 



2 

_ Y_ 
2 






r 

2 



dv = 



Jim' 



dt-y 






Eq. (4-8) becomes 



K 



-y e 



2 i 



e 

o Jim ' 



dt 



00 00 



K [ Pve 

^ J j — 7 ] ' dy + 

° ° Vl-P^ J2m 



1 1 
_ y_ 



K e 
' o 



°725r 



dt dy . 



where 



a a p 
K = -2LJ1: 

JUm 
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The first term in the above equation is zero; hence 
(4-8) becomes 



cr cr p^(l-p^)^ 
2ir 



+ K 



2 t 

r” - V f"^ 2 

I e — dt dy 



(4-9) 



Now considering the second term in (4-9) , 



» £ 
1 e- ^ 



-r -• 



_ Y- 






dt dy = J e ^ (p • 



5 + 



-r - 



o 



dt 



) dy 



(4-10) 

This result, due to the fact that y changes over positive 
numbers and p is between -.5 and zero, makes the argument 
of (•) positive. Hence (4-10) becomes 



, 


^ P 




2 


I 2 ^ . 


e 



2 

2 



dt dy = M 



Let 



t = V J2 

y = ^/2^? 



then 



M 



1 

~ 2 4 2 



-u f’-s 2 ^ 

— — J e ^ \f2 dv du 

o J2m o 



where 
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-s » 



2n 



M 



H -u /» zn 

ir’"" 



dv du 



Since the series is uniformly convergent, it fellows that 



M 



W VkW 

I p — n o 

1 I fr 1 I e \. .n 

‘2-J2*TZ\ 



^/2‘n 



(-P) 



2n+l 



o u 



n=0 (n/1-P ) 



n.i 



u 



^ 2n+l nl (2n+l) 



du 






or. 



n=0 



Iz£l 



2n+l 






pn 



2n+l 2n+l 



00 

f n -u 

J ^ ® 



du 



n: 



M 



2n+l 






n=0 



♦ 2n+l 



It was previously shown that p is always between -.5 and 0, 
Therefore, 

- 1 < - < 1 



Hence the series in the above expression is the expansion 
of Arc tan ( 7^ [Ref. 6] . 

Vi-p2 
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Thus Eq. (4-10) becomes 



1 / £ + _1 
2 V 2 ^ /- 



Arc tan 



n/2tt 



Ct%) 

s/l-p 



(4-11) 



Substituting (4-11) into (4-9) yields 
2 2 ^ 

<T <T p (1-p ) cap 



2TT 






Combining the results for 1, 2, and 3 into (4-7), and sub- 
stituting into (4-5) yields 



J J xy f(x,y)dy dx = 



2 — 

a cr (1-p ) 2 a a p 



- y X 



o o 



2TT 



+ + 



cr a p^(l-p^) cap 



2tt 



— C "7 J ~ ~ Arc tan (— 7===.)') 
y2rr ^ 2 V 2 y^_p2 ^ 



Thus, 



^ (j (j CTcTjO cr^ — A N 

(IxyI) = 4(^-^(lV)^ + ^ Arcta,n ( 7 =%)) - 

V 1 — P^ 



cap 

X y 



2 a a 



f^Cil-p^)^-p Arctan ("p^) ) (4-12) 

7l-p2 
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I I 

/ ^ 

But or^ and or^ are both equal to cr ^ 1 + * Hence, 



a 



E ( ! xyI ) = 



l+g ^ .2^h 



(d-p^) 



- p Arctan (—r=^=^)j (4-13) 

V 1— p 



Therefore, 



j n-2 2 20 ! 

V(A) = — cr - 



M. 



t-2 t-1 



(iV) »r(l+^)2 



a^+ 20(2 ^ i3i+:'E(I^_.e^_.l) 

i=0 j=d+l 



■3 

(4-14) 



where E ( | e^_ ^e^_^ I ) is given by (4-13). 

A 

It is clear that to arrive at an expression for V(A) 
in the form of (4-14) many tedious and difficult calcula- 
tions must be made. In fact, in order to arrive at an 

A 

explicit value for V(A), one must deal with an arctangent 

2 

term. However, because or is common to each term, the 
ratio V(A)/cr can be defined in terms of the known con- 
stants, and an approximate niimerical answer can be found 

•a 2 

for the ratio V(A)/cr . The computer program in Appendix 
B computes the first two terms in (4-13) exactly, and the 
third term approximately. Briefly, the computation of the 
third term was accomplished as follows: 

Eq. (4-13) is a bounded number for a given a and 
It cannot be made greater than some finite number d, what- 
ever be the value of p . (Note that p is the only variable 
in (4-13) that depends on i and j, and it is bounded.) 



36 



Therefore, in the double summation of (4-14) , for some 
finite value of i and j, higher order terms will get smaller 
than a specified number m which, in this case, was chosen 
as 10 The program was written so that the higher terms 

smaller than 10 were not included in the result, and 
were subsequently dropped from the summation. The upper 
bound of the error introduced by omitting the higher terms 
from the summation can be found as follows: 

E(1xy 1) is bounded from above by 

. = 4 ^ ^ _2_ 

ff(l+/9) 2 l+j9 

Since the summation is infinite. 

Upper Bound of Error = 



^ 10 -“ 



I I 

i=0 j=i+l 



40 :^ 10 ~^^ $ 



10 -“ 

(i+$> 



First, this error is additive and hence is added to 
the result. Second, whatever be the value of P , it is 
always much smaller than one. Hence, the results contained 
in Table 1, obtained by the computer program for various 
values of oi were almost exact. 
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SMOOTHING 

CONSTANT 


VARIANCE OF ( A) /ct^ 


STANDARD 

DEVIATION OF (or)/cr 


0.10 


0.0204 


0.1745 


0.15 


0.0325 


0.2173 


0.20 


0.0461 


0.2553 


0.25 


0.0614 


0.2905 


0.30 


0.0785 


0.3239 


0.35 


0.0979 


0.3561 


0.40 


0.1196 


0.3876 


0.45 


0.1440 


0.4187 


0.50 


0.1716 


0.4495 


0.55 


0.2026 


0.4802 


0.80 


0.4267 


0.6341 



Table 1. Variance of — •^r and Standard Deviation of — . 
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V. COMPARISON OF EXPONENTIAL SMOOTHING 



AND MLE BY COMPUTER SIMULATION 



For the reason explained previously, simulation was 
utilized to explore the parts of the model where theoreti- 
cal work was very difficult, and was also used to compare 
the exponential smoothing with MLE. The reason that MLE 
was chosen as the alternative procedure was because it has 
the property of an asymptotically vanishing variance. As 
indicated in [Ref. s], the maximum likelihood estimator. 



A 



<7 , of O' is given by the formula 



o 



A 




i=l 



where x is the sample mean. The statistic o has a chi- 
square distribution [Ref. 3 ]; in fact. 




and 




n§) ^2 



From these formulas. 




Lim 



V(o) = 0 
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The computer program used for the simulation is included in 
Appendix A. The simulation procedure and ensuing results 
are as follows : 

One thousand normally-distributed random numbers, re- 
presenting the demand, were generated by a pseudorandom 
number generator. Various values of the parameters, II, cr, 
and l/) were used. For every triplet of the parameters, re- 
order levels, risk levels attained, estimates of CT and 
sample mean, and sample standard deviation of c and of the 
estimate of cr by MAD were computed. The value of 0l was 0.2 
for all cases. Table 2 presents the results achieved by 
the simulation. As Table 2 indicates, the theoretical 
values of the standard deviations of both estimates of cr 
were in complete agreement with the sample standard devia- 
tions of these same estimates. The variance of or” is higher 

A 

than that of O’ in every case. The last two columns give 
the mean square error (MSE) of o and o. The MSE is given 

A 

by the variance plus the square of the bias. Although o 
is biased for 1000 observations, its bias can be neglected. 
Then the MSE of o and o are determined by their variances. 
The comparison of the MSE columns are quite striking; the 
MSE for o reaches values as high as 100 times that for O. 

By noting the similarity between the theoretical values and 
the simulation results, the effect on reorder levels was 
analyzed as follows: 

Computation of reorder levels was made by utilizing 
Eq. (3-3) . First, the theoretical values of the reorder 
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MSE 

C7 

\ 


.0004 

.0004 

.0016 

.0036 

! .0081 

! .0121 

( 

.0484 

.1236 

.4489 

.7921 

1.232 

1.795 

2.433 

3.920 

4.040 

4.272 


MSE 

a 


.065 

.065 

.260 

.585 

1.040 

1.625 

6.517 

26.01 

58.52 

104.04 

162.56 

234.09 

318.62 

416.16 

526.70 

651.78 


THEORETICAL 
SD OF $ 

i 

\ ■■ 


.02 

.02 

.04 

.06 

.03 

.11 

.22 

.44 

.67 

.89 

1.11 

1.34 

1.56 

1.98 

2.01 

2.23 


SAMPLE 
SD OF a 


.05 

.05 

.07 

.20 

.13 

.26 

.43 

.87 

1.51 

1.61 

2.42 

4.43 
5.01 
4.54 
4.85 
4.03 


THEORETICAL 
SD OF “a 


.255 

.225 

.510 

.765 

1.020 

1.275 

2.553 

5.100 

7.650 

10.200 

12.750 

15.300 

17.850 

20.400 

22.950 

25.530 


SAMPLE 
SD OFCT 


.26 

.26 

.50 

.75 

1.05 

1.36 
2.63 
5.81 

7.36 
10.40 
11.85 

15.43 
18.71 

19.44 
22.23 
25.10 


BIAS 

a 


-.04 

-.01 

-.07 

-.22 

.22 

.09 

-.01 

-.23 

-.67 

-.02 

-3.48 

-2.37 

-1.99 

-4.63 

-2.52 

-4.23 


BIAS 

a 


-.01 
.01 
-.05 
-.13 
.28 
.06 
.85 
.06 
.28 
1.23 
-3.09 
.20 
-1.58 
-2.4 3 
1.28 
-.80 


AVERAGE 

a 


.96 

.99 

1.93 

2.78 

4.22 

5.09 

9.99 

12.77 
29.33 
39.98 
46.52 
57.63 
68.01 
75.27 
87.48 

95.77 


AVERAGE , 
'a 


.99 

1.01 

1.95 

2.87 

4.28 

5.06 

10.35 

20.06 

30.28 
41.23 
46.91 
60.20 
67.42 
77.57 

91.28 
) 99.20 


PARAMETER i 
PAIRS 

(M, <y ) 


10,1 

10, 

20,2 

30.3 

40.4 

50.5 
100,10 
200,20 
300,30 
400,40 
500,50 
600,60 
700,70 
800,80 
900,90 

1000, IOC 



I 
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Table 2. Sample Characteristics 




levels were computed for every fixed parameter triplet )i, 
a, j|) . Then the reorder levels, found by smoothing and 
MLE, were compared with the theoretical values. This com- 
parison is summarized in Table 3. The value in Table 4 
is the percentage (rounded off to the nearest integer) of 
1000 reorder levels found to be within one unit of the 
theoretical reorder level. The first column under each of 
the headings P^q/ ^nd lists the results for smooth- 

ed estimates; the second columrj under each heading lists 
MLE. The quantity represents the number of units about 

the theoretical reorder level within which 50 per cent of 
the computed reorder levels were found; P^^ is a similar 
computation for 95 per cent of the reorder levels. Table 
3 clearly illustrates the superiority of MLE compared to 
exponential smoothing. For the case where smoothing was 
used to estimate cr, reorder levels fluctuated around the 
theoretical reorder level, causing either an excessive 
purchase, or imposing more risk than desired. The magni- 
tudes of these fluctuations are clearly evident from the 
table , 

Table 4 indicates the risk levels attained by both 
methods at the 1000th observation. For every predetermined 
parameter pair fi and cr, five different risk levels were 
selected. For every triplet, the theoretical reorder 
levels were compared with the computed reorder levels of 
both procedures. The actual risk levels attained were then 
computed. 
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PARAMETER 

TRIPLETS 


Qi 


^50 


^95 


40,4, .01 


27 


95 


2 


0 


6 


0 


40,4, .11 


41 


94 


2 


0 


4 


0 


40,4, .50 


49 


92 


1 


0 


3 


1 


50,5, .01 


22 


96 


3 


0 


8 


0 


50,5, .11 


32 


98 


2 


0 


5 


0 


50,5, .50 


46 


99 


1 


0 


3 


0 


100.10, .01 


12 


90 


5 


0 


14 


1 


100,10, .11 


18 


94 


4 


0 


11 


0 


100,10, .50 


26 


94 


2 


0 


6 


0 


200,20, .01 


5 


71 


10 


1 


29 


6 


200,20, .11 


9 


82 


6 


1 


18 


3 


200,20, .50 


10 


84 


5 


1 


14 


2 


300, 30, .01 


4 


81 


12 


1 


37 


6 


300,30, .11 


7 


74 


8 


1 


24 


3 


300,30, .50 


8 


51 


7 


1 


18 


3 


400,40, .01 


2 


15 


20 


3 


53 


10 


400,40, .11 


4 


32 


13 


2 


37 


6 


400,40, .50 


5 


75 


10 


1 


27 


3 


500,50, .01 


3 


1 


24 


9 


58 


11 


500,50, .11 


3 


4 


15 


4 


38 


6 


500,50, .50 


4 


73 


11 


1 


29 


4 


600,60, .01 


1 


12 


27 


3 


77 


17 


600,60, .11 


3 


10 


18 


3 


51 


5 


600 , 60 , . 50 


3 


1 


14 


4 


42 


12 


700,70, .01 


1 


33 


32 


3 


93 


16 


700,70, .11 


3 


35 


20 


2 


63 


10 


700,70, .50 


4 


16 


15 


2 


45 


15 


800,80, .01 


1 


0 


37 


12 


98 


27 


800,80, .11 


2 


0 


25 


8 


64 


18 


800,80, .50 


3 


32 


18 


2 


54 


9 



Table 3. Percentiles for Reorder 
Levels . 
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PARAMETER 

TRIPLETS 


Qi 


o 

in 


P 

95 


900,90, .01 


1 


30 


42 


2 


125 


12 


900,90, .11 


2 


37 


29 


2 


83 


12 


1000,100, .01 


1 


25 


45 


3 


126 


12 


1000,100, .11 


1 


8 


29 


4 


87 


7 


1000,100, .50 


1 


6 


25 


2 


68 


13 



Table 3 . (Continued) - Percentiles for Reorder 

Levels . 
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As Table 4 indicates, the risk levels attained by MLE 
are found to be very close to the desired risk levels. 
However, the risk levels of the smoothing method were 
found to fluctuate over a wide range. 

These results, gained by the simulation, indicated 
that the MLE technique is to be preferred over the tech- 
nique of exponential smoothing. However, due to the lack 
of theoretical analysis concerning the distribution of 
these results are not considered conclusive. However, this 
analysis does indicate that whenever a constant mean de- 
mand can be justified, the MLE technique is to be prefer- 
red to the exponential smoothing technique. 
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The first column under each heading lists 
exponential smoothing procedure. 



APPENDIX A 



SI MUL&TI ON 



C IN th= F'^LLOWING program two MF-^HQDS of F.STIM4TION OF 

C MFAN ANr> STANDARD OEVIATION 0«= A STATinNAPY ’^ROCFSS 

r APF CONTPASTFD BY USING SIMULATION t^chNIOUFS 

C method^ CRE: 

0 I.THF exponentially SMOQTHEO data, ts THE FSTTMiTF OF 

C Mcan and the exponentially SMOOTHEO ABSOIUTE VALUES 

C OF FPROPS AS THE ESTIMATE OF STamo/m?d DEVIATION TIMES 

C 4 KNOWN CONSTANT 

C 2 . maximum likelihood ESTIMATES 

'■ ALSO IN THE PROGRAM THE EFFECTS OF ROTH FSTIMATIQIvj 

C PROCEOUPFS ON THF RISK AND OPOPOEP LEVELS ARC <^imulA- 

C TED B Y S.‘MpLE SI 7E ICOO 



DIMENSION fK(6),T(6),Xl(6)TXP(6),XM(6)XN(6),XK(6I 
IXK PC0,5) , X2I 10G0,5) , G1 ( 10 JS ) »G2 ( 1 0 05 ) , 7 ( 10 ) » TX ( 6 ) 
2TYI6) ,PISKX(6) ,RI$KY(6) 



C RISK LEVELP 



PKI 1)=C.01 
OKI 2)=C.C5 
PK( 3)=C. 11 
OKI 4)=c.25 
5) = C. fC 

READI5.131) (T(I),I=1,E) 
131 format! 5F4. 2) 

S=10.0 
XJ=C.C 
Y=URN( C) 
no e21 MT=1,16 
XJ=XJ+1.C 
DO 1 J=l,5 
DO 1 1=1,800 
XII I, J I =0.0 
1 X2fI,J)=C.C 
7L=S 



C smoothing constant alpha is 3.2 



BETA^O.P 

WRITE(6,07E) S,STDEX 

Q7C cqpma T( • 1 'MEAM=* ,nX, •P.D'=V=' ,F6. 1 ) 

DELTA =C.7Fe*SORT(2.C/(l.L+BETA) )*STOi=X 



C COMOllTA TION CF THE THEORETICAL REORDER LEVELS 



no 11 1=1,5 
11 7( I )= 7L + STDEX*T(I ) 

WRTTE(6,15C) (7(1) ,T=1,5) 

15C FORMA T( //TIG, 'THEORETICAL REOROEP LEVELS • , 5 F 10 . 3/ / ) 
V = 1.G 
SUM=7L 
TnT=o.C 
R 1=0.0 
B2=0. 0 
B3=C.O 
R4=C.O 



C generation of normally distibuteo numbers which have 
C mcan s and ptandard deviation STDEX 
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DO 21 1=1,1000 

CAIL NORMAL ( C. 0 , STDF X , X > 

Y=7L+X 



c method 1 



DELTA=0.2»ABS( S-Y) + C. 8 *DELTA 

S=C.2*Y+C. S*S 

RDEV=1. 18PTDELT6 

V=V+1.C 

S(]M = 5UM + Y 



C method 2 



X8AR=SUM/v 

TQT=TOT+( XBAO-Y) **2 

STDEV=S0PT( T0T/( V-l.C) ) 

OK n = SDEV 

G2( I ) = STDEV 

DO 20 J=l,5 

XL( J ) = S + SOE V*T( J» 

XP( J)=XRAR+STDE V*T( J) 
XK( J)=Xl ( J)-Z( J) 

TX( J )=( XL( J)-ZL ) /STDEX 



C RISKS ATTAINED BY BOTH METHODS 



RTSKXI J)=C.5-0.5TERF(TX(J) *U.7C7) 

TY( J)=( XP( JI-ZL) /STDFX 
R ISKY( J)=C. 5-0. 5*EPP ( TY( J) *0.7C7) 

XM( .|)=XP( J)-Z( J) 

M1=XK( J) 

N2=XM( J) 

IF(Nl.LT.O) N1=-M1 
IP(N2.LT.O) N2=-N2 
XKNl + 1, J) = X1(N1+1 ,J)+1.0 

20 X2»N2+1,J)=X2(N2+1 ,J)+1.0 
NS=I 

IF( NS .LE.SRP) GO TC 21 

WRTTF(6,1C1) ( XL( M) ,M = 1 ,5) , ( XP( V) , M=i , 5 ) 

ICl FORMAT! /T13 , 5F7. 1 , lOX ,5F7. 1 ) 

WRITF(6,147) (RISKX(M) ,M=1 ,5) ,(RISKY(M> ,M=1,5) 
147 FORMAT! /T2, 'ACTUAL P T SK • ,5 F7 . 3 , lOX , 5 F7 . 3 ) 

21 CONTINUE 



C COMPUTATION CF THE OBSERVED FREQUENCIES WITHIN THE I 

C NEIGHBORHOOD OF ACTUAL REORDER LEVEL WHERE 1=1, 2, 3,.. 



DO 26 J=l,5 
WRITF!6,105) PK!J) 

1G5 FORMAT! //T45 , *P I SK = • , F5 . 2 ) 

A 1=0.0 
A2=C.O 
1 = 1 

25 IF!A1 .GE. l.OI GO TO 24 
A 1=X1! I , J) /ICOO.O+Al 
24 IF! A 2 .GE. 1.0) GO TO 35 
A2=X2! I ,J) /1C00.0+A2 
35 WRITE!6,112) I,A1,A2 

112 FORMAT! Tio, tPERCENTILE BETWEEN thEO. REORDER LEVEL 
1 +AND- • , I4,2F10.4) 

1 = 1 + 1 
XA=A 1+A2 
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IF( X4 .LF. 1.996J GOTO 25 
26 CONTINUF 



C sample and STANDAPD deviation OF BOTH ESTIMATES 

C OF POPULATION STANDARD DEVIATION 



DO 27 I=l,lCOO 
B1=G1( I )+Bl 

27 B2=82+G2< I ) 

B1=B1/1C0C.0 
B2=B2/1CCC.0 
DO 2B 1=1,1000 
B3=B3+(G1( I )-Bl »*»2 

28 8A=BA + (B2-G2n n**2 
B3=SORT(B3/1COO.O) 

B4=SORT(p.4/ic0O.0) 

WRITE! 6t200) 

200 FORMAT! //TA5, ’SAMPLE MEAN AND S. DEV OF ESTIMATED S. 
IDEV’ ) 

WRITE!6»2Cn B1,83,B2»RA 

2C1 '=OPMA T! //TAB, M MAD) • ,5X, ’MFAN=’ ,F7 . 2 , • S DEV= • , «=7.2, 15X, 
1 •! S2) ’ ,fX, 'MEAN=’ ,F7.2, 'SDEV=' ,F7.2) 

S=10C.C 

IF(XJ-5.C) 721,721 ,722 

721 S=XJ*1C.0 
GO TO 82 1 

722 S=S«! XJ-5.0) 

821 continue 

STOP 

END 



SUBROUTINE NORMAL !EX,STDX,X» 

SUM=O.C 

DO 5 1=1,26 

R=URM! 1) 

5 SUM=SUM+P 
$U«=SUM/26.0 

X=STDX*17.61*! SUM-0.5)+EX 

RETURN 

END 
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APPENDIX B 



NIJVEOICAI COMPUT/tTION OF VA°MNCF CF A 



IMPLICIT R^AL^fl (/-H,C-Z) 
^FTi= . o 

n(i 20 K=i,ii 



C SMOPTHTNO CCNFTANT 



ALPHA = 1.-B5TA 

A FUM = o.C 
WRTTCI 6 , 110 ) 

110 FORMA T( //TIC ,* 'SMOOTH! NG CONSTANT = • , F^ . 3 ) 



C APCRC Xiwtf TION of THF INFINITE F!JM 



no 1 M= 1,303 

1=M-1 

J=M 

T-4.0/( 3.1'a 1*( l.O+BFTAM 
P = Al oHA/( l. + BET/ ) 
n=RPTA««( I+J) 

C=RE'^A«’M >)-I ) 

A=C /BE^A 
3 R=A«BE'^A 

FHO = ( 1. ♦BET^. ) /■- . ♦( 1. + F*( 1. u-A’X^Z + I 1 .+B ) **2) + ( 1.- ALPHA* 
lA )**2- 1. 

7N=nSQR''( 1 ,-RHn**2 ) 

7C=-rhO/7.N 

RL)M = T*( 7N-RHO*nATAN( 7P) ) *0 
A Fijy=A £t'M+ 'l)v 
I<=( SlJM-1. CP-10) 1,1,2 
2 A=A*BETA 
n= 0 #BETA 
GO TO 3 
1 CONTiNi.'E 

WRITF(6 102) SUM 

1C? FQPMfl T( TIC, ' LAST SUM= ' ,E20. 12 ) 

X=ALDHA/( 1 .4-RFTA ) **2*2. 2 82 /3. l‘-l 
Z=P .C /3. 141*RFTA /( 1.0+RETA) **2 



C VARIANCE QS A 



VAP=X-Z + 2.0*( 1. 0-RETA ) **2*aS.iJM 



c ST/NDARn n-VIA''ION OR a 



Rnpv=0 sop'll 3.1-^1*(1. + B'^TA) *VAR) /2.C 
WR!TE(6,1CI) VAR,SOFV 

ICl PORMA TITIO, 'VARI ANCE MA 0 • , FI 5 . 7 , 5 X , * S AND ART OEVIA 
ITTON SFIE.?) 

BE'^A = PFTA-.05 
IRIK. SO. 1C) BETA=.2 
20 CONTINUE 
RTOP 
END 
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